********* Home-specific cost-benefit analyses

********************************
**** average monthly retail electricity prices
clear all
import delimited "$maindir/Data/electricity_prices", rowrange(11) delimiter(comma)

** correctly identifying month and year of observation
forval j = 1/12 {
	local year = 2001
	local start = `j'+3
	forval i = `start'(12)215{
		rename v`i' m`j'y`year'
		local year = `year' + 1
	}
}

drop v1 v2 v3
drop *y2001 *y2002 *y2003 *y2004 *y2005 *y2006 *y2007 *y2008 *y2017 *y2018


** calculate average for each month
forval j = 1/12 {
	egen avg_m`j' = rowmean(m`j'*)
	sum avg_m`j'
	sca elec_price_m`j' = r(mean)
}


********************************
**** average monthly retail natural gas prices
clear
import excel "$maindir/Data/retail_gas.xls", sheet("Data 1") cellrange(A3:B386) firstrow

rename IllinoisPriceofNaturalGasDe price
gen year = year(Date)
gen month = month(Date)
drop Date
drop if price==.

keep if year>=2009 & year<=2016
forval j = 1/12 {
	sum price if month==`j'
	sca gas_price_m`j' = r(mean)
}


********************************
**** average monthly citygate natural gas prices
clear
import excel "$maindir/Data/citygate_gas.xls", sheet("Data 1") cellrange(A3:B386) firstrow

rename NaturalGasCitygatePriceinIl price
gen year = year(Date)
gen month = month(Date)
drop Date
drop if price==.

keep if year>=2009 & year<=2016
forval j = 1/12 {
	sum price if month==`j'
	sca gas_citygate_m`j' = r(mean)
}


***************************
**** incorporating social cost of carbon

/**** social marginal costs of electricity
from by Borenstein and Bushnell (2018): https://www.nber.org/papers/w24756
dataset provided by the authors */

** retail electricity prices by provider
clear 
use "$maindir/Data/retail_final.dta"

keep if eia_id_d == 56697
keep year res_sales avgcharge

sum avgcharge [aweight = res_sales] if year==2014, detail
sca avg_price2014 = `r(mean)'
sum avgcharge [aweight = res_sales] if year==2015, detail
sca avg_price2015 = `r(mean)'
sum avgcharge [aweight = res_sales] if year==2016, detail
sca avg_price2016 = `r(mean)'

** correcting for social costs
clear 
use "$maindir/Data/BBCorrection_UtilityData.dta"

sort date hour
gen month = month(date)
gen month_year = mofd(date)
format month_year %tm
order month_year, after(year)

gen retail_prices = avg_price2014 if year==2014
replace retail_prices = avg_price2015 if year==2015
replace retail_prices = avg_price2016 if year==2016

* monthly social marginal costs
bysort month_year : egen smc_monthly = mean(smc)

* monthly sales
bysort month_year : egen res_sales_monthly = total(res_sales)

* drop duplicates
keep month month_year retail_prices smc_monthly res_sales_monthly varcharge
duplicates drop

gen pricediff = smc_monthly - retail_prices

* monthly price differences
forval j = 1/12 {
	sum pricediff [aweight = res_sales_monthly] if month==`j'
	sca elec_pricediff_m`j' = `r(mean)'
}



******************************************* loading the cost-benefit data
clear
use "$maindir/Results/Model_Outputs/cba_data.dta"

**** single observation per month per home
keep Household end_month weight_boots* monthsave_post* ///
	monthuse_counterpost* prismrest rsquare_pre rsquare_post hdd_pre hdd_post
drop if monthsave_post==.
duplicates drop

egen nboots = rownonmiss(weight_boots*)
drop if nboots==0
duplicates drop

gen aux = 1
by Household : egen nmonths = total(aux)
order Household end_month nmonths, first
sort Household end_month 



********************************************************************************
**** monetizing the benefits

************************************
* price escalation based on "Energy Price Indices and Discount Factors..."
* https://nvlpubs.nist.gov/nistpubs/ir/2018/NIST.IR.85-3273-33.pdf
* Table Ca-2
* assuming no more escalation after 30 years

** escalation for electricity
local counter = 0
foreach x in 1.02 1.07 1.09 1.09 1.09 1.09 1.10 1.10 1.09 1.09 1.09 1.09 1.10 1.10 1.10 ///
	1.10 1.09 1.09 1.09 1.09 1.09 1.08 1.08 1.08 1.07 1.07 1.07 1.07 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 {
		local counter = `counter'+1
		sca avg_elec_doll_mmbtu = 0
		forval j = 1/12 {
			sca elec_cents_mmbtu = (elec_price_m`j' + 10*elec_pricediff_m`j') / 0.003412 /* monthly average from 2009-2016, converting to MMBtu */
			sca elec_doll_mmbtu = elec_cents_mmbtu/100
			sca avg_elec_doll_mmbtu = avg_elec_doll_mmbtu + elec_doll_mmbtu
			sca elec_doll_mmbtu`counter'_m`j' = `x'*elec_doll_mmbtu
		}
}
sca avg_elec_doll_mmbtu = avg_elec_doll_mmbtu/12
di "Average electricity price = `=round(scalar(avg_elec_doll_mmbtu),0.01)'"
	
** escalation for natural gas
/*
price of carbon: $40 per ton CO2
emissions factor from ton CO2 per Mcf: 0.0549
damages from carbon: $2.196 per Mcf of natural gas


********************************************************************************
************* damages from other pollutants
https://www.pnas.org/content/pnas/suppl/2019/04/03/1816102116.DCSupplemental/pnas.1816102116.sapp.pdf
Table S3 - Fuel combustion: natural gas – residential 

* emissions factors
https://www3.epa.gov/ttnchie1/conference/ei12/area/haneke.pdf
https://19january2017snapshot.epa.gov/sites/production/files/2015-08/documents/eiip_areasourcesnh3.pdf - Table III

damages from SO2: $25.65 per ton
emissions factor SO2 lb/MMBtu: 0.6/1020
emissions factor SO2 ton/MMBtu: (0.6/1020)/2204.62 = 0.000000266882
damages from SO2: $0.0000068 per MMBtu

damages from NOx: $17.67 per ton
emissions factor NOx lb/MMBtu: 94/1020
emissions factor NOx ton/MMBtu: (94/1020)/2204.62 = 0.0000418
damages from NOx: $0.0007 per MMBtu

damages from PM2.5: $358.92 per ton
emissions factor PM2.5 lb/MMBtu: 1.9/1020
emissions factor PM2.5 ton/MMBtu: (1.9/1020)/2204.62 = 0.000000844
damages from PM2.5: $0.0003 per MMBtu

damages from Amonia: $105.43 per ton
emissions factor Amonia lb/Mcf: 20/1000
emissions factor Amonia ton/Mcf: 20/1000/2204.62 = 0.000095
damages from Amonia: $0.01 per Mcf

other pollutants do not cause significant damages, due to their low emissions factors
********************************************************************************

final price of carbon: monthly citygate price + 2.196 
*/
local counter = 0
foreach x in 1.03 1.06 1.08 1.09 1.12 1.13 1.14 1.15 1.16 1.16 1.17 1.17 1.17 1.18 1.18 ///
	1.19 1.19 1.20 1.21 1.22 1.22 1.23 1.23 1.24 1.25 1.25 1.26 1.27 1.28 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 {
		local counter = `counter'+1
		sca avg_gas_doll_mmbtu = 0
		forval j = 1/12 {
			sca gas_doll_mmbtu = (gas_citygate_m`j' + 2.196) / 1.037   /* monthly average from 2009-2016, converting to MMBtu */
			sca avg_gas_doll_mmbtu = avg_gas_doll_mmbtu + gas_doll_mmbtu
			sca gas_doll_mmbtu`counter'_m`j' = `x'*gas_doll_mmbtu
		}
}
sca avg_gas_doll_mmbtu = avg_gas_doll_mmbtu/12
di "Average gas price = `=round(scalar(avg_gas_doll_mmbtu),0.01)'"

**** monetizing the benefits in each year
sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */

forval i = 1/80 {
	gen elecbenefits_year`i' = .
	forval j = 1/12 {
		replace elecbenefits_year`i' = (-1)*monthsave_post*elec_contrib*elec_doll_mmbtu`i'_m`j' if end_month==`j'
	}
}

forval i = 1/80 {
	gen gasbenefits_year`i' = .
	forval j = 1/12 {
		replace gasbenefits_year`i' = (-1)*monthsave_post*(1-elec_contrib)*gas_doll_mmbtu`i'_m`j' if end_month==`j'
	}
}

forval i = 1/80 {
	gen totbenefits_year`i' = elecbenefits_year`i' + gasbenefits_year`i'
	drop elecbenefits_year`i' gasbenefits_year`i'
}


**** bootstrap monetized benefits
egen savings_SD = rowsd(monthsave_post*)
gen min95_monthsave = monthsave_post + 1.96*savings_SD /* note that min and max signs invert because negative values denote enery savings */
gen max95_monthsave = monthsave_post - 1.96*savings_SD

sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */

foreach x in min95 max95 {
	forval i = 1/80 {
		gen elecbenefits_year`i' = .
		forval j = 1/12 {
			replace elecbenefits_year`i' = (-1)*`x'_monthsave*elec_contrib*elec_doll_mmbtu`i'_m`j' if end_month==`j'
		}
	}

	forval i = 1/80 {
		gen gasbenefits_year`i' = .
		forval j = 1/12 {
			replace gasbenefits_year`i' = (-1)*`x'_monthsave*(1-elec_contrib)*gas_doll_mmbtu`i'_m`j' if end_month==`j'
		}
	}

	forval i = 1/80 {
		gen `x'benefits_year`i' = elecbenefits_year`i' + gasbenefits_year`i'
		drop elecbenefits_year`i' gasbenefits_year`i'
	}
}


*** finally aggregate to household level
bysort Household : egen realized_savings_mmbtu = total(monthsave_post)
forval j = 1/200 {
	bysort Household : egen realized_savings_mmbtu`j' = total(monthsave_post`j')
}

bysort Household : egen counterfactual_mmbtu = total(monthuse_counterpost)
forval j = 1/200 {
	bysort Household : egen counterfactual_mmbtu`j' = total(monthuse_counterpost`j')
}

forval i = 1/80{
	bysort Household : egen hhbenefits_year`i' = total(totbenefits_year`i')
	bysort Household : egen min95_hhbenefits_year`i' = total(min95benefits_year`i')
	bysort Household : egen max95_hhbenefits_year`i' = total(max95benefits_year`i')
}


drop end_month monthsave_post* monthuse_counterpost* totbenefits_year* min95benefits_year* max95benefits_year* savings_SD min95_monthsave max95_monthsave
duplicates drop



****** CE for each home
* cost share for each upgrade category (excluding H&S)
merge m:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(tot_act* nonHS_cost Real_nonHS_cost) keep(3) nogen
drop if Real_nonHS_cost==.

foreach x in Attic AirCon AirSeal Baseload Door WtHtr Foundation General Window Furnace WallIns {
gen `x'_costshare = tot_act`x'/nonHS_cost
}
egen sharecheck = rowtotal(*costshare)


** cost-weighted average of lifetime of retrofits for each home
* based on WeatherWorks documentation: "Table SIR1 - Retrofit Service Lives"
gen avg_lifetime20 = 20*AirSeal_costshare + 25*WallIns_costshare + 25*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 25*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime20 = round(avg_lifetime20)
gen avg_lifetime_months20 = avg_lifetime20*12
replace avg_lifetime_months20 = round(avg_lifetime_months20)

* assuming all measures have a 10-year lifespan
gen avg_lifetime10 = 10*AirSeal_costshare + 10*WallIns_costshare + 10*Attic_costshare + ///
	10*Window_costshare + 10*Door_costshare + 10*Baseload_costshare + 10*AirCon_costshare + ///
	10*Furnace_costshare + 10*Foundation_costshare + 10*General_costshare + 10*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime10 = round(avg_lifetime10)
gen avg_lifetime_months10 = avg_lifetime10*12
replace avg_lifetime_months10 = round(avg_lifetime_months10)

* insulation measures with 50 years lifespan
gen avg_lifetime30 = 20*AirSeal_costshare + 50*WallIns_costshare + 50*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 50*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime30 = round(avg_lifetime30)
gen avg_lifetime_months30 = avg_lifetime30*12
replace avg_lifetime_months30 = round(avg_lifetime_months30)

* insulation measures with 80 years lifespan
gen avg_lifetime40 = 20*AirSeal_costshare + 80*WallIns_costshare + 80*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 80*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime40 = round(avg_lifetime40)
gen avg_lifetime_months40 = avg_lifetime40*12
replace avg_lifetime_months40 = round(avg_lifetime_months40)



******** Program cost-effectiveness by amount spent
sca elec_contrib = 0.17

**** 30 years lifetime
sca rate = 0.03
gen year30_benefits = 0
gen min95_year30_benefits = 0
gen max95_year30_benefits = 0
forval i = 1/50 {
	gen year30_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace year30_benefits = year30_benefits + year30_benefits`i'
	drop year30_benefits`i'
	
	gen year30_benefits`i' = min95_hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace min95_year30_benefits = min95_year30_benefits + year30_benefits`i'
	drop year30_benefits`i'

	gen year30_benefits`i' = max95_hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace max95_year30_benefits = max95_year30_benefits + year30_benefits`i'
	drop year30_benefits`i'
}
gen year30_npv = year30_benefits - Real_nonHS_cost
gen min95_year30_npv = min95_year30_benefits - Real_nonHS_cost
gen max95_year30_npv = max95_year30_benefits - Real_nonHS_cost


**** no discount rate
sca rate = 0
gen disc0_benefits = 0
forval i = 1/50 {
gen disc0_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace disc0_benefits`i' = 0 if round_lifetime30<`i' & disc0_benefits`i'!=.
replace disc0_benefits = disc0_benefits + disc0_benefits`i'
drop disc0_benefits`i'
}
gen disc0_npv = disc0_benefits - Real_nonHS_cost


**** 6% discount rate
sca rate = 0.06
gen disc6_benefits = 0
forval i = 1/50 {
gen disc6_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace disc6_benefits`i' = 0 if round_lifetime30<`i' & disc6_benefits`i'!=.
replace disc6_benefits = disc6_benefits + disc6_benefits`i'
drop disc6_benefits`i'
}
gen disc6_npv = disc6_benefits - Real_nonHS_cost


**** 40 years lifetime
sca rate = 0.03
gen year40_benefits = 0
forval i = 1/80 {
gen year40_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year40_benefits`i' = 0 if round_lifetime40<`i' & year40_benefits`i'!=.
replace year40_benefits = year40_benefits + year40_benefits`i'
drop year40_benefits`i'
}
gen year40_npv = year40_benefits - Real_nonHS_cost


**** 20 years lifetime
sca rate = 0.03
gen year20_benefits = 0
forval i = 1/25 {
gen year20_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year20_benefits`i' = 0 if round_lifetime20<`i' & year20_benefits`i'!=.
replace year20_benefits = year20_benefits + year20_benefits`i'
drop year20_benefits`i'	
}
gen year20_npv = year20_benefits - Real_nonHS_cost


**** 10 years lifetime
sca rate = 0.03
gen year10_benefits = 0
forval i = 1/10 {
gen year10_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year10_benefits = year10_benefits + year10_benefits`i'
drop year10_benefits`i'	
}
gen year10_npv = year10_benefits - Real_nonHS_cost


**** export to be used in other analyses
keep Household realized_savings_mmbtu* counterfactual_mmbtu* weight_boots* Real_nonHS_cost year30_benefits-year10_npv nmonths prismrest rsquare_pre rsquare_post hdd_pre hdd_post
duplicates drop
label data ""

save "$maindir/Results/Model_Outputs/cba_scc.dta", replace
export delimited using "$maindir/Results/Model_Outputs/cba_scc.csv", replace




********************************************************************************
********************************************************************************
********************** retail energy prices

******************************************* loading the cost-benefit data
clear
use "$maindir/Results/Model_Outputs/cba_data.dta"


**** single observation per month per home
keep Household end_month weight_boots* monthsave_post* ///
	monthuse_counterpost* prismrest rsquare_pre rsquare_post hdd_pre hdd_post
drop if monthsave_post==.
duplicates drop

egen nboots = rownonmiss(weight_boots*)
drop if nboots==0
duplicates drop

gen aux = 1
by Household : egen nmonths = total(aux)
order Household end_month nmonths, first
sort Household end_month 



********************************************************************************
**** monetizing the benefits

************************************
* price escalation based on "Energy Price Indices and Discount Factors..."
* https://nvlpubs.nist.gov/nistpubs/ir/2018/NIST.IR.85-3273-33.pdf
* Table Ca-2
* assuming no more escalation after 30 years

** escalation for electricity
local counter = 0
foreach x in 1.02 1.07 1.09 1.09 1.09 1.09 1.10 1.10 1.09 1.09 1.09 1.09 1.10 1.10 1.10 ///
	1.10 1.09 1.09 1.09 1.09 1.09 1.08 1.08 1.08 1.07 1.07 1.07 1.07 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 ///
	1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 1.06 {
		local counter = `counter'+1
		sca avg_elec_doll_mmbtu = 0
		forval j = 1/12 {
			sca elec_cents_mmbtu = (elec_price_m`j') / 0.003412 /* monthly average from 2009-2016, converting to MMBtu */
			sca elec_doll_mmbtu = elec_cents_mmbtu/100
			sca avg_elec_doll_mmbtu = avg_elec_doll_mmbtu + elec_doll_mmbtu
			sca elec_doll_mmbtu`counter'_m`j' = `x'*elec_doll_mmbtu
		}
}
sca avg_elec_doll_mmbtu = avg_elec_doll_mmbtu/12
di "Average electricity price = `=round(scalar(avg_elec_doll_mmbtu),0.01)'"


** escalation for natural gas
local counter = 0
foreach x in 1.03 1.06 1.08 1.09 1.12 1.13 1.14 1.15 1.16 1.16 1.17 1.17 1.17 1.18 1.18 ///
	1.19 1.19 1.20 1.21 1.22 1.22 1.23 1.23 1.24 1.25 1.25 1.26 1.27 1.28 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 ///
	1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 1.29 {
		local counter = `counter'+1
		sca avg_gas_doll_mmbtu = 0
		forval j = 1/12 {
			sca gas_doll_mmbtu = (gas_price_m`j') / 1.037   /* monthly average from 2009-2016, converting to MMBtu */
			sca avg_gas_doll_mmbtu = avg_gas_doll_mmbtu + gas_doll_mmbtu
			sca gas_doll_mmbtu`counter'_m`j' = `x'*gas_doll_mmbtu
		}
}
sca avg_gas_doll_mmbtu = avg_gas_doll_mmbtu/12
di "Average gas price = `=round(scalar(avg_gas_doll_mmbtu),0.01)'"

**** monetizing the benefits in each year
sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */

forval i = 1/80 {
	gen elecbenefits_year`i' = .
	forval j = 1/12 {
		replace elecbenefits_year`i' = (-1)*monthsave_post*elec_contrib*elec_doll_mmbtu`i'_m`j' if end_month==`j'
	}
}

forval i = 1/80 {
	gen gasbenefits_year`i' = .
	forval j = 1/12 {
		replace gasbenefits_year`i' = (-1)*monthsave_post*(1-elec_contrib)*gas_doll_mmbtu`i'_m`j' if end_month==`j'
	}
}

forval i = 1/80 {
	gen totbenefits_year`i' = elecbenefits_year`i' + gasbenefits_year`i'
	drop elecbenefits_year`i' gasbenefits_year`i'
}


**** bootstrap monetized benefits
egen savings_SD = rowsd(monthsave_post*)
gen min95_monthsave = monthsave_post + 1.96*savings_SD /* note that min and max signs invert because negative values denote enery savings */
gen max95_monthsave = monthsave_post - 1.96*savings_SD

sca elec_contrib = 0.17 /* assuming 17% of savings are from electricity */

foreach x in min95 max95 {
	forval i = 1/80 {
		gen elecbenefits_year`i' = .
		forval j = 1/12 {
			replace elecbenefits_year`i' = (-1)*`x'_monthsave*elec_contrib*elec_doll_mmbtu`i'_m`j' if end_month==`j'
		}
	}

	forval i = 1/80 {
		gen gasbenefits_year`i' = .
		forval j = 1/12 {
			replace gasbenefits_year`i' = (-1)*`x'_monthsave*(1-elec_contrib)*gas_doll_mmbtu`i'_m`j' if end_month==`j'
		}
	}

	forval i = 1/80 {
		gen `x'benefits_year`i' = elecbenefits_year`i' + gasbenefits_year`i'
		drop elecbenefits_year`i' gasbenefits_year`i'
	}
}


*** finally aggregate to household level
bysort Household : egen realized_savings_mmbtu = total(monthsave_post)
forval j = 1/200 {
	bysort Household : egen realized_savings_mmbtu`j' = total(monthsave_post`j')
}

bysort Household : egen counterfactual_mmbtu = total(monthuse_counterpost)
forval j = 1/200 {
	bysort Household : egen counterfactual_mmbtu`j' = total(monthuse_counterpost`j')
}

forval i = 1/80{
	bysort Household : egen hhbenefits_year`i' = total(totbenefits_year`i')
	bysort Household : egen min95_hhbenefits_year`i' = total(min95benefits_year`i')
	bysort Household : egen max95_hhbenefits_year`i' = total(max95benefits_year`i')
}


drop end_month monthsave_post* monthuse_counterpost* totbenefits_year* min95benefits_year* max95benefits_year* savings_SD min95_monthsave max95_monthsave
duplicates drop


****** CE for each home
* cost share for each upgrade category (excluding H&S)
merge m:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(tot_act* nonHS_cost Real_nonHS_cost) keep(3) nogen
drop if Real_nonHS_cost==.

foreach x in Attic AirCon AirSeal Baseload Door WtHtr Foundation General Window Furnace WallIns {
gen `x'_costshare = tot_act`x'/nonHS_cost
}
egen sharecheck = rowtotal(*costshare)


** cost-weighted average of lifetime of retrofits for each home
* based on WeatherWorks documentation: "Table SIR1 - Retrofit Service Lives"
gen avg_lifetime20 = 20*AirSeal_costshare + 25*WallIns_costshare + 25*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 25*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime20 = round(avg_lifetime20)
gen avg_lifetime_months20 = avg_lifetime20*12
replace avg_lifetime_months20 = round(avg_lifetime_months20)

* assuming all measures have a 10-year lifespan
gen avg_lifetime10 = 10*AirSeal_costshare + 10*WallIns_costshare + 10*Attic_costshare + ///
	10*Window_costshare + 10*Door_costshare + 10*Baseload_costshare + 10*AirCon_costshare + ///
	10*Furnace_costshare + 10*Foundation_costshare + 10*General_costshare + 10*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime10 = round(avg_lifetime10)
gen avg_lifetime_months10 = avg_lifetime10*12
replace avg_lifetime_months10 = round(avg_lifetime_months10)

* insulation measures with 50 years lifespan
gen avg_lifetime30 = 20*AirSeal_costshare + 50*WallIns_costshare + 50*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 50*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime30 = round(avg_lifetime30)
gen avg_lifetime_months30 = avg_lifetime30*12
replace avg_lifetime_months30 = round(avg_lifetime_months30)

* insulation measures with 80 years lifespan
gen avg_lifetime40 = 20*AirSeal_costshare + 80*WallIns_costshare + 80*Attic_costshare + ///
	15*Window_costshare + 15*Door_costshare + 15*Baseload_costshare + 15*AirCon_costshare + ///
	20*Furnace_costshare + 80*Foundation_costshare + 10*General_costshare + 15*WtHtr_costshare ///
	if sharecheck>0.98
gen round_lifetime40 = round(avg_lifetime40)
gen avg_lifetime_months40 = avg_lifetime40*12
replace avg_lifetime_months40 = round(avg_lifetime_months40)



******** Program cost-effectiveness by amount spent
sca elec_contrib = 0.17

**** 30 years lifetime
sca rate = 0.03
gen year30_benefits = 0
gen min95_year30_benefits = 0
gen max95_year30_benefits = 0
forval i = 1/50 {
	gen year30_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace year30_benefits = year30_benefits + year30_benefits`i'
	drop year30_benefits`i'
	
	gen year30_benefits`i' = min95_hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace min95_year30_benefits = min95_year30_benefits + year30_benefits`i'
	drop year30_benefits`i'

	gen year30_benefits`i' = max95_hhbenefits_year`i'/((1+rate)^`i')
	replace year30_benefits`i' = 0 if round_lifetime30<`i' & year30_benefits`i'!=.		
	replace max95_year30_benefits = max95_year30_benefits + year30_benefits`i'
	drop year30_benefits`i'
}
gen year30_npv = year30_benefits - Real_nonHS_cost
gen min95_year30_npv = min95_year30_benefits - Real_nonHS_cost
gen max95_year30_npv = max95_year30_benefits - Real_nonHS_cost


**** no discount rate
sca rate = 0
gen disc0_benefits = 0
forval i = 1/50 {
gen disc0_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace disc0_benefits`i' = 0 if round_lifetime30<`i' & disc0_benefits`i'!=.
replace disc0_benefits = disc0_benefits + disc0_benefits`i'
drop disc0_benefits`i'
}
gen disc0_npv = disc0_benefits - Real_nonHS_cost


**** 6% discount rate
sca rate = 0.06
gen disc6_benefits = 0
forval i = 1/50 {
gen disc6_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace disc6_benefits`i' = 0 if round_lifetime30<`i' & disc6_benefits`i'!=.
replace disc6_benefits = disc6_benefits + disc6_benefits`i'
drop disc6_benefits`i'
}
gen disc6_npv = disc6_benefits - Real_nonHS_cost


**** 40 years lifetime
sca rate = 0.03
gen year40_benefits = 0
forval i = 1/80 {
gen year40_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year40_benefits`i' = 0 if round_lifetime40<`i' & year40_benefits`i'!=.
replace year40_benefits = year40_benefits + year40_benefits`i'
drop year40_benefits`i'
}
gen year40_npv = year40_benefits - Real_nonHS_cost


**** 20 years lifetime
sca rate = 0.03
gen year20_benefits = 0
forval i = 1/25 {
gen year20_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year20_benefits`i' = 0 if round_lifetime20<`i' & year20_benefits`i'!=.
replace year20_benefits = year20_benefits + year20_benefits`i'
drop year20_benefits`i'	
}
gen year20_npv = year20_benefits - Real_nonHS_cost


**** 10 years lifetime
sca rate = 0.03
gen year10_benefits = 0
forval i = 1/10 {
gen year10_benefits`i' = hhbenefits_year`i'/((1+rate)^`i')
replace year10_benefits = year10_benefits + year10_benefits`i'
drop year10_benefits`i'	
}
gen year10_npv = year10_benefits - Real_nonHS_cost


**** export to be used in other analyses
keep Household realized_savings_mmbtu* counterfactual_mmbtu* weight_boots* Real_nonHS_cost year30_benefits-year10_npv nmonths prismrest rsquare_pre rsquare_post hdd_pre hdd_post
duplicates drop
label data ""

save "$maindir/Results/Model_Outputs/cba_retail.dta", replace
export delimited using "$maindir/Results/Model_Outputs/cba_retail.csv", replace



